function [G0, G1, C, PSI, PI, pos_c, pos_pc, pos_mc, pos_yk, pos_p, pos_ps, pos_zk]  = Gensys_matrix(n, frequency, omeg_c, Omeg, params)

% number of sectors
n_sectors = n;

% sectoral Calvo
FPA = frequency;
alpha = (ones(size(FPA)) - FPA);

% sectoral weights in consumption
Omega_C = omeg_c;

%sectoral weights 
Omega = Omeg;


%% Parameters

% Parameters values
phi_pi = params.phi_pi;          % param in Taylor rule
phi_c = params.phi_c;            % param in Taylor rule
% new parameters for monetary policy rule
rho_i = params.rho_i;
phi_gc = params.phi_gc;
rho_a = params.rho_a;            % persistence agg tech shocks
rho_mu = params.rho_mu;          % persistence monetary shock
rho_k = params.rho_k;            % persistence idiosync tech shock
beta = params.beta;              % discount factor
frisch = params.frisch;          % inverse of Frisch elasticity
sigma = params.sigma;            % CRRA pameter
eta = params.eta;                % elast of subst across sectors
theta = params.theta;            % elast of subst within sectors
delta = params.delta;            % delta in production function

% compute the compound parameters

% psi in production function 
psi = delta*(theta - 1)/theta;   


invk_k = zeros(1,n_sectors);
for i = 1:n_sectors
    invk_k(i) = alpha(i)/((1 - alpha(i))*(1 - beta*alpha(i)));
end


% sectoral participation in production
n_k = (1-psi).*inv(eye(n_sectors)-psi.*Omega)*Omega_C;

% sectoral outdegrees
outdegree = 1/psi.* (n_k - (1-psi)*Omega_C);

% predefine variables
pos_ps = zeros(n_sectors,1);
pos_p = zeros(n_sectors,1);
pos_Ep = zeros(n_sectors,1);
pos_mc = zeros(n_sectors,1);
pos_w = zeros(n_sectors,1);
pos_yk = zeros(n_sectors,1);
pos_zk = zeros(n_sectors,1);
pos_lk = zeros(n_sectors,1);
pos_ak = zeros(n_sectors,1);

% Position indices for the variables
pos_c = 1;
pos_Ec = 2;
pos_pc = 3;
pos_Epc = 4;
pos_mu = 5;
pos_y = 6;
pos_a = 7;
pos_i = 8;
maxind = 8;


for i = 1:n_sectors
    
    % superscript k prices
    pos_ps(i) = maxind + i;
    % subscript k variables
    pos_p(i) = maxind + n_sectors + i;
    pos_Ep(i) = maxind + 2*n_sectors + i;
    pos_mc(i) = maxind + 3*n_sectors + i;
    pos_w(i) = maxind + 4*n_sectors + i;
    pos_yk(i) = maxind + 5*n_sectors + i;
    pos_zk(i) = maxind + 6*n_sectors + i;
    pos_lk(i) = maxind + 7*n_sectors + i;
    pos_ak(i) = maxind + 8*n_sectors + i;
end
maxind2 = 9;


%% Equations
% cast into coefficient matrices for Gensys:
% G0*y(t) = G1*y(t-1) + C + PSI*z(t) + PI*eta(t)

% matrices G0, G1, C, PSI and PI. include nominal interest rate now (the
% +1 relative to previous code)
G0 = zeros(maxind+maxind2*n_sectors, maxind+maxind2*n_sectors);
G1 = zeros(maxind+maxind2*n_sectors, maxind+maxind2*n_sectors);
C  = zeros(maxind+maxind2*n_sectors, 1);
PSI= zeros(maxind+maxind2*n_sectors, 2+n_sectors); 
PI = zeros(maxind+maxind2*n_sectors, 2+n_sectors);

id1 = 0;
id2 = 0;

% IS:
% ct = Ect+1 - 1/sigma * (it-Epct+1 +pct)
id1 = id1 + 1;
G0(id1,pos_c) = 1;
G0(id1,pos_Ec) = -1;
G0(id1,pos_i) = 1/sigma;
G0(id1,pos_Epc) = -1/sigma;
G0(id1,pos_pc) = 1/sigma;

% Taylor rule:
id1 = id1 + 1;
% it = rhoi it-1+(1-rhoi)(phi_c ct + phipi pit + phigc (ct-ct-1));
G0(id1, pos_i)=1;
G1(id1, pos_i)=rho_i;
G0(id1, pos_c)=-(1-rho_i)*(phi_c+phi_gc);
G1(id1, pos_c)=-(1-rho_i)*phi_gc;
G0(id1, pos_pc)=-(1-rho_i)*phi_pi;
G1(id1, pos_pc)=-(1-rho_i)*phi_pi;
G0(id1,pos_mu) = -1; 


% price quasi PCs
for k = 1:n_sectors
    G0(id1+k,pos_p(k)) = -(1/alpha(k) + alpha(k)*beta)*invk_k(k);
    G0(id1+k,pos_Ep(k)) = beta*invk_k(k);
    G0(id1+k,pos_mc(k)) = 1;
    G1(id1+k,pos_p(k)) = -invk_k(k);
end
id2 = id2 + 1;


% marginal costs
for k = 1:n_sectors
    G0(id1+ id2*n_sectors + k, pos_mc(k)) = 1;
    G0(id1+ id2*n_sectors + k, pos_w(k)) = -(1-delta);
    G0(id1+ id2*n_sectors + k, pos_ps(k)) = -delta;
    G0(id1+ id2*n_sectors + k, pos_ak(k)) = 1;
    G0(id1+ id2*n_sectors + k, pos_a) = 1;
end
id2 = id2 + 1;

% wages
for k = 1:n_sectors
    G0(id1+ id2*n_sectors + k, pos_w(k)) = 1;
    G0(id1+ id2*n_sectors + k, pos_pc) = -1;    
    G0(id1+ id2*n_sectors + k, pos_c) = -sigma;
    G0(id1+ id2*n_sectors + k, pos_lk(k)) = -frisch;
end
id2 = id2 + 1;

for k = 1:n_sectors
    G0(id1+ id2*n_sectors + k, pos_w(k)) = 1;
    G0(id1+ id2*n_sectors + k, pos_ps(k)) = -1;    
    G0(id1+ id2*n_sectors + k, pos_zk(k)) = -1;
    G0(id1+ id2*n_sectors + k, pos_lk(k)) = 1;
end
id2 = id2 + 1;

% sectoral output demand side 
for k = 1:n_sectors
    G0(id1+ id2*n_sectors + k, pos_yk(k)) = n_k(k);
    G0(id1+ id2*n_sectors + k, pos_c) = -(1-psi)*Omega_C(k);
    G0(id1+ id2*n_sectors + k, pos_pc) = -(1-psi)*Omega_C(k)*eta;
    G0(id1+ id2*n_sectors + k, pos_p(k)) = psi*outdegree(k)*eta + (1-psi)*Omega_C(k)*eta;
    for k2 = 1:n_sectors
        G0(id1+ id2*n_sectors + k, pos_zk(k2)) = -psi*n_k(k2)*Omega(k,k2);
        G0(id1+ id2*n_sectors + k, pos_ps(k2)) = -psi*n_k(k2)*Omega(k,k2)*eta;
    end
end
id2 = id2 + 1;

% sectoral output - supply side
for k = 1:n_sectors
    G0(id1+ id2*n_sectors + k, pos_yk(k)) = 1;
    G0(id1+ id2*n_sectors + k, pos_ak(k)) = -1;
    G0(id1+ id2*n_sectors + k, pos_a) = -1;
    G0(id1+ id2*n_sectors + k, pos_lk(k)) = -(1-delta);
    G0(id1+ id2*n_sectors + k, pos_zk(k)) = -delta;
end
id2 = id2 + 1;  



% defining aggregate prices pc, ptilde, pk
%aggregate output 
id1 = id1 + 1; 
G0(id1+ id2*n_sectors, pos_y) = 1;
for k=1:n_sectors
    G0(id1+ id2*n_sectors, pos_yk(k)) = -n_k(k);
end
% y
id1 = id1 + 1;
G0(id1+ id2*n_sectors, pos_pc) = 1;
for k=1:n_sectors
    G0(id1+ id2*n_sectors, pos_p(k)) = -Omega_C(k);
end
% ps
for k=1:n_sectors
    G0(id1+ id2*n_sectors + k, pos_ps(k)) = 1;
    for k2 = 1:n_sectors
        G0(id1+ id2*n_sectors + k, pos_p(k2)) = -Omega(k2,k);
    end
end
id2 = id2 + 1;

% monetary shock process [12]
id1 = id1 + 1;
G0(id1+ id2*n_sectors, pos_mu) = 1;
G1(id1+ id2*n_sectors, pos_mu) = rho_mu;
PSI(id1+ id2*n_sectors, 1) = 1;


% technology shock process
id1 = id1 + 1;
G0(id1+ id2*n_sectors, pos_a) = 1;
G1(id1+ id2*n_sectors, pos_a) = rho_a;
PSI(id1+ id2*n_sectors, 2) = 1;

for k = 1:round(n_sectors)
    G0(id1+ id2*n_sectors+ k, pos_ak(k)) = 1;
    G1(id1+ id2*n_sectors+ k, pos_ak(k)) = rho_k;
    PSI(id1+ id2*n_sectors+ k, 2+k) = 1;
end
id2 = id2 + 1;

id1 = id1 + 1;
G0(id1+ id2*n_sectors, pos_c) = 1;
G1(id1+ id2*n_sectors, pos_Ec) = 1;
PI(id1+ id2*n_sectors, 1) = 1;

id1 = id1 + 1;
G0(id1+ id2*n_sectors, pos_pc) = 1;
G1(id1+ id2*n_sectors, pos_Epc) = 1;
PI(id1+ id2*n_sectors, 2) = 1;

for k = 1:n_sectors
    G0(id1+ id2*n_sectors+ k, pos_p(k)) = 1;
    G1(id1+ id2*n_sectors+ k, pos_Ep(k)) = 1;
    PI(id1+ id2*n_sectors+ k, 2+k) = 1;
end
id2 = id2+1;

G0(isinf(G0)) = 0;  
G0(isnan(G0)) = 0;        
G1(isinf(G1)) = 0;  
G1(isnan(G1)) = 0; 